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ABSTRACT 

We consider local, stratified, numerical models of isothermal accretion disks. 
The novel feature of our treatment is that radial extent and azimuthal extent 
Ly satisfy H -C Lx,Ly <^ R, where H is the scale height and R is the local 
radius. This enables us to probe mesoscale structure in stratified thin disks. We 
evolve the model at several resolutions, sizes, and initial magnetic field strengths. 
Consistent with earlier work, we find that the saturated, turbulent state consists 
of a weakly magnetized disk midplane coupled to a strongly magnetized corona, 
with a transition at \z\ ~ 2H. The saturated a ~ 0.01 — 0.02. A two-point 
correlation function analysis reveals that the central AH of the disk is dominated 
by small scale turbulence that is statistically similar to unstratified disk models, 
while the coronal magnetic fields are correlated on scales ~ lOH. Nevertheless 
angular momentum transport through the corona is small. A study of magnetic 
field loops in the corona reveals few open field lines and predominantly toroidal 
loops with a characteristic distance between footpoints that is ~ H. Finally 
we find quasi-periodic oscillations with characteristic timescale ~ 30f2~^ in the 
magnetic field energy density. These oscillations are correlated with oscillations 
in the mean azimuthal field; we present a phenomenological, alpha-dynamo model 
that captures most aspects of the oscillations. 

Subject headings: accretion, accretion disks, magnetic fields, corona, magnetohy- 
drodynamics 



1. Introduction 

The physics of angular momentun i transport is at th e core of accretion disk stud- 



i es. C l assical viscous thin d i sk the ories ( IShakura fc SunyaevI ( Il973l ) ; iLynden-Bell fc Pringle 



( ll974j ): iNovikove fc Thornd ( 119731 )) assume the existence of a local turbulent viscous stress. 
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thus provide a simple local parameterization, i.e., "anomalous viscosity" a, for disk momen- 
tum transport and dissip ation. Since the early 90's, the magnetorotational instability (MRI, 
Balbus fc Hawleyl ( 1199 ll )) has been regarded as the best candid ate to drive accretion disk 
turbulence, although gravitational torque or magnetic winds of a lBlandford fc Paynd (Il982l ) 
type can also enhance angular momentum transport. 

Classical thin disk theories are vertically integrated and azimuthally averaged, there- 
fore essentially one dimensional. Currently, disk vertical structure can only be obtained 
from numerical simulations where turbulence is established from first-principle instabili- 
ties such as the MRI. Glo b al disk simul a tions a re just start i ng to investigate thin disks 
faevnolds fc FabianI J2008h : Ishafee et all J2008h : iNoble et al.l J2009I )). but they are com- 
putationally expensive and not yet able to fully resolve turbulent structures in the disk. 
Shearing box simulations, on the other hand, can concentrate resolution on disk dynamics 
at scales of order the disk scale height H = c^/i^, therefore are more suitable to study ac- 
cretion flows in detail. Pa s t studies of s hearing box simulations with vertical gravity (e.g ., 



Brandenburg et a. 



(119951): 



Stone et al. ( 1996 ): Miller k Stond (l2000h : iHirose et al.l (120061 ): 



Blaes et al.l (120071 ): ISuzuki fc Inutsukal ( l2009l )) have revealed a rich set of structures and 
dynamics in stratified disks. However, all these stratified shearing box simulations were 
done with a box of limited radial extent ~ H, therefore t hey w ere not able to explore 
any structure on scales larger than H. Rece ntlv. iDavis et al. (120101) have studied stratified 
shearing box of radial extent = AH, and iJohansen et al.l (|2009[ ) have adopted models of 
box size up to ~ lOH in their zonal flow studies. However both these studies are limited 
to the small veritical extent (~ ±2H) and physically unrealistic periodic vertical boundary 
conditions. In this paper we study the dynamics and structure in isothermal stratified disks 
using large shearing box with domain sizes > lOH in all directions. 

We still do not know whether a magnetized turbulent disk is well modeled as a steady- 
state, locally dissipated disk model. It is possible, for example, that structures (gas and/or 
fields) develop at a scale large compared to H, and that these structures could be associ- 
ated with nonlocal energy or angular momentum transport. Large scale structures might 
also develop in the magnetic field i n the f orm o f dynamo. The disk might also be secu- 
larly unstable (see the overview by iPiranI (Il978l )). that could cause the disk to break up 
into rings. It is well known that a Navier-S tokes viscosity model for d isk turbulence leads, 
for so me opacity regimes, to both viscous (ILightman fc Eardleyl 1 19741 ) and thermal iPiran 
(119781 ) instability, although it is now believed that thermal instability can be removed by 



delays imposed through finite relaxation time effects in MRI-driven turbulence (IHirose et al. 
J2009h ). 



^From an observational point of view, the level of fluctuations (inhomogeneity) at dif- 
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ferent locations in disks and how these different locations communicate with each other have 



(Davis et al. 


2005; 


Blaes et al. 


2006) 



these models observational diagnostics require integrating over the disk surface, so radially 
extended structure in the disk model may change the disk spectrum. Our disk model is 
isothermal (we do not solve an energy equation) and is therefore not capable of investigating 
dissipation and radiation. It is possible that larger fluctuations would appear in physi- 
cally richer rnodels w here thermody namics and radiative effects are taken into account (e.g., 



Turner et al.l (120031 ): iTurneii (120041 )). It would then be interesting in the future for spectral 



modelers to consider disk models with larger radial domains. 

A shearing box larger than H is also esse ntial to catch the field structure and dynamics 
in the accretion disk magnetic coronae (ADC: lTout fc Pringld (119961 ): also see a discussion in 
Uzdensky fc GoodmanI (120081 ) ). where the field has a characteristic curvature I ~ fa/^ > H, 
and Va is the characteristic Alfven speed in the region. 

Recently, it has also been pointed out that a large box size may be important to study 
the saturatio n properties of the MRI-dri ven turbulence, either on the ground of resolving par- 
asitic mo des (IPessah fc Goodman! |2009| ). or in a phenomenological model of an MRI driven 
dynamo ( Vishnia3 20091 ). Saturation mechanisms in stratified disk may be fundamentally 
different from those in unstratified disks. Recent numerical experiments on unstratified disks 
suggest that: (a) with a zero- net flux, the saturation is depende nt on the microscopic Prandtl 
number Pr^ in the disk, at least at the low Reyn olds number (IFromang &: Papaloizoull2007l : 
Lesur fc LongarettilboO?! : Isimon fc HawlevljgOO^ : (b) with a net (toroidal or ve rtical) flux, 
the saturation increases with resolution (IHawley et al.l Il995l : iGuan et al.l l2009l ). Stratified 
disk models, which are closer to real disks, may well maintain a net (most likely, toroidal; 
see a discussion in iGuan fc Gammid (120091 )) field in the disk region because of the mag- 
netic buoyancy induced by stratification. Therefore we expect a saturation in stratified disk 
models to differ from unstratified models. 

It is worth enumerating the assumptions we adopt in this work: (1) we use an isothermal 
equation of state (EOS) in our models; (2) the vertical support comes from the gas and 
magnetic pressure rather than the radiation pressure; (3) there is no explicit viscosity or 
resistivity; (4) our initial conditions consist of a uniform toroidal field in a region near the 
disk midplane; (5) we use outflow boundary conditions for the vertical boundaries. 

The paper is organized as follows. In §2 we give a description of the local model and 
summarize our numerical algorithm. In §3 we present a fiducial model and analyze its 
structure in the saturated state. In §4 we describe how this structure depends on model 
parameters. In §5 we give a report on quasiperiodic oscillations ( "butterfly diagrams" ) and 
present a phenomenological model to describe them that is based on a mean-field dynamo 
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model; §5 contains a summary of our results. 



2. Local Model and Numerical Methods 

The local model for disks can be obtained by expanding the equations of motion around a 
circular-orbiting coordinate origin, with (r, 0, z) = (ro, Qot + <Po, 0) in cylindrical coordinates, 
assuming that the peculiar velocities are comparable to the sound speed and that the sound 
speed is small compared to the orbital velocity. The local Cartesian coordinates are then 
obtained from cylindrical coordinates via (x, y, z) = {r — Tq, ro[0 — Qot — 4>o]: z)- this work 
we assume the disk sits in a Keplerian (1/r) potential. We also use an isothermal {p = cj.p, 
where Cg is constant) EOS. 

For an ideal MHD disk, the equation of motion in the local model is 

^ + v.Vv + cl^ + ^-^-^-^ + 2nxv-3n'xx + n'zz = 0. (1) 
dt p Svrp Anp ^ ' 

/^From left to right, the last three terms in Eqn ([1]) represent the Coriolis force, tidal forces 
and vertical gravitational acceleration in the local frame respectively. The orbital velocity 
in the local model is ^ 

Vorh=-l^'^V- (2) 

This velocity, along with a vertical density profile p{z) = po exp[— fi^z^/(2c^)] and zero 
magnetic field, is a steady-state solution to Eqn(IT]). po is the midplane density. In this 
work, we nondimensionalize the local model by choosing po = 1, = 1, and Cg = 1; the 
usual disk scale height H is therefore H = Cg/fl = I. The initial surface density is therefore 
J pdz = \/2ttPo. 

The local model is realized numerically using the "shearing box" boundary conditions 
(e.g. iHawley et al.lll995l ). which isolates a rectangular region in the disk. The azimuthal (y) 
boundaries are periodic; the radial (x) boundaries are "shearing periodic" ; they connect the 
radial boundaries in a time-dependent way that enforces the mean shear. The vertical (z) 
boundaries use a form of outflow boundary conditions: all variables in ghost zones (including 
the z velocity and momentum on vertical boundaries because of the staggered mesh) are 
copied from the last active zone in the computational domain, with the additional constraint 
that no inflow is allowed. For stratified disk models the outflow boundary condition is better 
motivated than periodic boundary conditions, although it is more difficult to implement. 



What constraint do these boundary conditions place on the field evolution? Integrating 
the induction equation over the computational domain yields, after application of Stokes 
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theorem, 



L^LyL^dt{B^) = dt d^xB^ = dx d[s] ■ {v x B 



(3) 



where the second integral is taken on a circuit round the box boundaries at fixed x. It is 
evident that the EMF integrated over a fine on the top boundary will not cancel that on 
the bottom boundary for outfiow boundary conditions, and so {B^) is not conserved. A 
similar argument implies that (By) is not conserved either. dt{Bz) is proportional to a line 
integral around the box at constant z, where the quasi-periodic radial and periodic azimuthal 
boundary conditions do cause cancellation, so (B^) is constant (numerically: constant to 
within accumulated roundoff error). 



In the preceding paragraph we adopted the notation ( ) for a volume average: 



[f) = — J dxdydzf. 



We will also use 

for a plane average, and 

for a time average. 



[/] 



LxLy 



dxdyf 



f I dtf 



(4) 
(5) 
(6) 



field (jjohnson et al 



20081: 



Our models are evolved using ZEUS (IStone fc Normanlll992l ) with "orbital advection" 



dMassetlboOol: ICammielboOli I Johnson Sz Gammid l2005l aka FARGO; see) for the magnetic 



Fromang fc Stondl2009l ). ZEUS is an operator-split, finite differ- 



ence scheme on a staggered mesh that uses a Von Neumann- Richtmyer artificial viscosity to 
capture shocks (this is a nonlinear bulk viscosity that does not produce significant angular 
momentum transport in our models), and the Method of Characteristics-Constrained Trans- 
port (MOC-CT) scheme to evolve the magnetic field and preserve the V-S = constraint to 
machine precision. The orbital advection is implemented on top of ZEUS . It decomposes the 
velocity field into a mean shear part with orbital velocity Vorb = —(l^x[y] and a fiuctuating 
part 6v] V = 6v + Vorb- Advection for the mean fiow can is done using interpolation (which 
is always stable), so that the Courant limit on the timestep depends only on 6v and not v^rb- 
Shearing boxes with > H, where the shear fiow is supersonic, can then be evolved more 
accurately, and with a larger timestep. 

We have also implemented an additional procedure to make th e numerical diffusion mor e 
nearly translation invariant in the plane of the disk. As discussed in lGuan &: Gammid (120091 ). 
the entire box is shifted by a few grid points in the radial direction at t = 2nLy/ {"iVLL^), 



- 6 - 



n = 1,2,3, . . .); at these instants the box is exactly periodic. After the shift we execute a 
divergence cleaning procedure to remove the monopoles that are created by joining the radial 
boundaries together in the middle of the computational domain. This procedure carries little 
computational cost. 

The timestep in large stratified disk simulations is limited through the Courant condition 
by the Alfven speed va = B/^/47^p at large \z\/H, where the density is orders of magnitude 
smaller than at 2; = 0. To prevent the simulation from being brought to a halt by low density 
zones (and to avoid other numerical artifacts associated with small p), we impose a density 
floor pmin = 10~^po- This density floor is ~ 1 — 2 orders of magnitude smaller than the 
averaged minimum density in the saturated state. We have tested a smaller density floor 
Pmin = 10"''po and found the choice of the density floor does not affect our results. 



3. Large Stratified Disk Simulations 

3.1. Fiducial Model 

All models start from a hydrodynamical equilibrium, with p{z) = po exp{—z'^ /[2H'^]) . 
We introduce a uniform toroidal field Bq = B^y at 1^1 < 2iJ; Bq is chosen so that at the disk 
midplane the initial plasma parameter /3o = SttPq/Bq = 25 (the sharp vertical variation in By 
at 1 2; I = 2H makes the disk initially unstable to magnetic Rayleigh- Taylor instability, but this 
structure is quickly wiped out by MRI driven turbulence). Each component of the velocity 
is perturbed in each zone, with 6vi uniformly distributed in [— 0.01, O.Oljc^,. The models are 
evolved long enough (> 150 orbits ~ 900i7~^) to reach a saturated, i.e., statistically steady, 
state. 

Our fiducial model has a domain size of {L^, Ly, L^) = (16, 20, 10)if and resolution 
384 X 256 X 128. This corresponds to a physical resolution of (24, 12.8, 12.8) zones per scale 
height. Snapshots of p and Eb = B'^/Sir at slices with constant x,y and z in the saturated 
state are shown in Figure [H 

Turbulence is confined to the region \z\ < 2H. Within this region magnetic field fluc- 
tuations are contained on a scale / -C H, with a structure in the shape of narrow filaments 
that are extended by the azimuthal shea r. This turbulent field structure resembles that 



observed in unstratified disk simulations (IGuan fc Gammid 120091 ). Density fluctuations on 



a scale ~ H due to sound waves are also evident in the x — z plane density snapshots. At 
\z\ > 2H the MRI is suppressed. Eb decreases sharply, but not as rapidly as p. 



The disk vertical structure is shown in Figure |2], which shows [p], [Eb], [/3], Maxwell 
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(magnetic) stress [M^y] = [—B^By]! At: and Reynolds stress [Rxy] = [pVxSvy]. These profiles 
are obtained from a time average over the last lOOvr/fi. The most striking feature in these 
profiles is the "turbulent disk surface" at \z\ ~ 2H defined by [Eb]{z) and [Mxy]{z). Inside 
this surface both are independent of z; outside both exhibit an approximately exponential 
dependence on z. As illustrated in the vertical profile of f3, as \z\ increases, magnetic energy 
density drops slower than density; above \z\ ~ 2.5H (5 drops below unity. Therefore the 
region \z\ > 2H is magnetically dominated and this leads to the suppression of the MRI. 
Prom now on we will simply refer to the magnetically dominated upper region with (3 < 1 
as "corona", and the turbulent \z\ < 2H region as "disk." 

Fits to the disk structure give 



,2 



and 



' 0.036po exp (— 44^^ ), otherwise. 
O.OUpocl if \z\ < 2.55H; 



[EB]iz) 



0.012 exp (—^Ll^)poc^, otherwise. 



0.64H 

In the saturated state, [p]{z) is different from the initial density profile poexp[—z'^/{2H'^)] 
due to the magnetic buoyancy effects and mass loss through the z boundaries. Inside the 
disk, a nearly Gaussian density profile indicates that this region is still mainly supported by 
gas pressure. 

How can we understand the vertical magnetic structure of the disk? A uniforna l y mag- 



netiz ed atmosphere is subject to interchange and Parker type modes flNewcomblll961t iParker 



19661 ). The more dangerous of these is Parker, whose stability condition is 



"^■^ gas 

where Pgas is the gas pressure, q = fl'^ z is the gravitational acceleration, and 7 is the adiabatic 



index (here, 7 = 1) (jNewcomblll96ll ). Por a disk in hydrostatic equilibrium 

pg; (10) 



dz 

together these conditions imply 

0. (11) 



dPmag dEs 



dz dz 



Marginally stable stratification therefore corresponds to constant [-E_b], as is found at \z\ < 
2H . This suggests that (1) magnetic buoyancy is driving the disk toward a marginally stable 
state; (2) magnetic buoyancy is crucial in controlling the vertical magnetic structure in the 
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bulk of the disk. If this is correct, it follows that Eb{z) in the disk could be different in 
nonisothermal models. In particular, marginal stability requires 

^'^\B'^-,p(l.^jBl^.'lhl\. (12) 



8tt dz \7 dz dz 

Thus an isentropic disk has dB^/dz = 0, while a stably stratified, nonradiative disk (in the 
Schwarzschild sense) can support dB'^/dz < 0. In a radiative disk the instability criterion is 
modified (disks heated by internal dissipation of turbulence rather than external irradiation 
tend to have strong radiative diffusion, or Peclet numbers of order ~ 50), because radial 
radiative diffusion tends to wipe out temperature perturbations for the most unstable modes 
with high radial wavenumbers. 

Figure [3] shows the evolution of magnetic energy density in the disk {Eb,a)i magnetic 
energy density in the corona {Eb,c) and (a) 

, < fW^yd^x , , 

J pCgd-^x 

where W^y = Rxy + M^y is the total shear stress Averaging the last 50 orbits, we found that 
H ~ 0.013, jE^/pocl ~ 0.012, jE^/pocl ~ 0.0043. 



3.2. Two-point correlation function 

One question motivating this study was whether thin disks exhibit mesoscale structure, 
i.e. structure on scales that are ^ H but ^ R. As is evident in Figure [H the characteristic 
scale of the magnetic energy density varies with \z\. Near \z\ = 0, turbulent structure 
resembles that observed in unstratified disk models: the field is confined in small structures 
with scale / ^ H. Away from \z\ = 0, / increases, reaching ~ i7 at |z| ~ 2.5H. 

The two-point correlation function ^ provides a quantitative measure of disk structure: 
^Biz) = [SB{x, y- z) ■ 6B{x + Ax,y + Ay; z)]. (14) 



Here S B = B-\B 



; for a detailed discussion of ^ and the corresponding correlation lengths 



Xi see (iGuan et al.ll2009l ). Figure H] shows ^b{z) in the {Ax, Ay) plane ai z = 0, z = 2.5H 
and z = A.bH. In these plots we have averaged 8 neighboring vertical zones to increase the 
signal-to-noise ratio. At the disk midplane, the correlation function has a narrow elliptical 
core of width A < a few H. As \z\ increase to 2.5H the core becomes larger, especially 
in the radiation direction, and low amplitude features develop on scales of ~ lOH. These 
low-amplitude, mesoscale features are new and are not seen in unstratified disk models. 
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3.3. Coronal loop structure 

Our disk models contain a "corona", where /3 < 1. It is not clear how accurately, or 
inaccurately, our code models this region because it contain s no explicit model for recon- 



nection (nor is any convincing model currently available; see lUzdensky &: Goodman! (120081 ) 
for a discussion of the difficulties of simulating force-free coronae). Still, it is interesting 
to characterize the field structure in existing simulations before asking how they might be 
changed by more sophisticated reconnection models. 

How can we understand coronal magnetic field structure? Most of the coronal field 
is anchored in the disk, so we begin by sampling field lines that rise through the surface 
z = 2.5H at a single instant. Using bilinear interpolation for the field, we trace field lines 
initiated from every cell on the z = 2.5H surface, until they either (a) come back to the 
z = 2.5H surface, or (b) leave the upper z surface, or (c) exceed maxmum integration step 
10^ indicating the formation of a closed loop. A snapshot of these field lines are show in 
Figure O Two features of the coronal field are obvious just from visual inspection: many of 
the field lines return to the disk after only a short sojourn in the corona, and the loops tend 
to have greater azimuthal than radial extent. 

A more quantitative appr oach is to calculate a coronal l oop distribution function, as in 



the phenomenological model of lUzdensky &: Goodman! (120081 ) (hereafter UG). The field lines 



should then be sampled according to the flux through each zone surface d^i = BzAdxdy. 
We also average over the last 50 orbits to improve the loop statistics. We find that ~ 96% 
of the field lines passing through the z = 2.5H surface return to the same surface, ~ 4% of 
field lines are open in the sense that the escape through the upper boundary of the box, and 
only ~ 0.1% of the field lines form closed loops inside the corona. We have also found that 
the small fraction of the open field lines is quite stable during the saturated state, ranging 
from ~ 2 — 5% at instaneous state, therefore it appears that the corona fields structure has 
reached a statistical steady state. 

We then use three variables to describe the geometry of close field lines (field loops) that 
return to the z = 2.5H surface: the loop foot point separation Ax in the x — y plane, the loop 
maximum height Azmax, and the loop orientation angle 6'foot = the angle between the foot 
separation vector and the y axis. We calculate the distributions functions d^/dAx, d^/dAy, 
d^/dAz^^^ and d^/dOfoot by following the trajectory of every field line that emerges at the 
center of each zone surface i, then weighting the result by the flux The final distribution 
function is normalized by |$|, the total absolute flux through the z = 2.5H plane. 

Figures |6] and [7] show d^/dAx, d^/dAy, d^/dAz^^,:^ and d^/dOfoot, averaged over the 
last 50 orbits. From d^/d9foot (Figure E]) it is evident that most of the field loops are 



- 10 - 



orientated in the azimuthal direction with 6 < a few degrees. This suggests that shear plays 
a significant role in determining the coronal field structures. Most loops also have maximum 
height Azmax < H (see d^/dOioot in Figure [7]). 

If there is no reconnection at all then magnetic energy injection from the underlying 
turbulent disk might cause the loop to grow in an unlimited way. This is not the case here: 
although we do not include dissipation explicitly, numerical reconnection due to truncation 
errors is present in our numerical scheme, as in all finite-difference MHD schemes. It is 
difficult, however, to quantify the numerical reconnection rate in our model directly. We 
therefore try to compare our numerical loop distributions to the predictions of the UG 
model. 

We fit a power-law to (i$ / rfAa; and d^ / dAy, 

where Co is a constant. Notice the shape for d^/dAx and d^/dAy are very similar. For 
Ay a fit between 3H < Ay < 20H gives k ^ —1.2. We can also calculate the general loop 
distribution function for AL = (Ax^ + Ay^)^/^. It almost overlaps with the d^/dAy curve 
because the loops are nearly toroidal, and a fit between 3H < Ay < 20H gives k ^ —1.2. 

In the model of UG, a slope of ~ —2 corresponds to the limit that reconnection is 
slow compared to the shear (the dimensionless reconnection parameter k, ~ 0.01 in UG), 
and a slope of k > —1.5 corresponds to the cases when the total magnetic energy of the 
corona is dominated by the largest loops {k < 0.002). The shallow k ~ —1.2 slope measured 
here then indicate our numerical models are probably in a slow Sweet-Parker reconnection 
regime. However, this comparison should not be taken too seriousljll] because our model is 
ideal MHD, and does not explicitly model reconnection. One serious concern is that the 
coronal reconnection could fall into a fast, coUisionless regime which is poorly understood, 
and not well modeled by our grid scale dissipation. 

Lastly we want to comment on several surface effects in our model. These include 
the yz component of magnetic stress tensor Myz = —ByBz/{4:7i), the vertical components 
of kinetic flux and Poynting flux, and the mass loss rate. Notice these quantities do not 
necessarily average to zero because of the outflow boundary conditions. We have found 
that (Myz) is nearly zero with temporal fluctuations of amplitude < 10~'^pqcI, much smaller 
than the dominant xy component [M^y). In the steady state, the vertical energy flux is 



^Although the reconnection rate in the corona might well determine the vertical magnetic energy profile 
[Eb,c{z)\, simply from a characteristic field curvature argument, where l{z) ~ [wq(z)]/J7 ^ Aa;(z). 
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dominated by the advective part of the Poynting flux, which is on the order of lO'^poc^. 
This vertical energy flux is only ~ 1% of the turbulent dissipation rate Q in the disk main 
body {Q ~ apoc^ ~ lO'^poC^), indicating a weak vertical energy flux. 

The disk loses mass through the upper and lower boundaries. In the last 50 orbits, the 
disk lost ~ 1.4% of its initial mass. The vertical mass loss rate is not negligibl^ in our 
Lz = lOH models because of the outflow boundary conditions. However, we have noted 
a trend in which the vertical mass loss decreases with increasing L^. For example, in our 
Lz = 12H model the disk lost ~ 0.64% of its initial mass during the last 50 orbits, giving a 
mass loss rate half that of the Lz = lOH model. We therefore expect a decreasing mass loss 
rate as Lz increases. 

It is also worth noting that in our models the mean vertical magnetic fields maintain 
{Bz) = because of the shearing box boundary conditions. We also found that the plane- 
averaged [Bz] ~ at all z including at the domain boundaries. The vertical field at the 
surfaces are turbulent with patches of opposite sign field penetrating the boundaries. How- 
ever the vertical field here are fluctuations with radial correlation length < a few H and 
amplitude ~ one order of magnitude smaller than that at the disk mid plane. We have not 
observed a steady magnetic wind and the observed mass loss is probably due to outflow 
boundaries. 

To summarize, in our models we see a weak wind launched from the disk surface. In a 
steady state both vertical energy and momentum flux are negligible. 



Here we give a brief discussion of the saturation dependence on model parameters, 
including: (1) resolution; (2) L^; (3) L^; (4) initial field strength in terms of plasma parameter 
Pq. When exploring parameter space we vary only one parameter at a time unless stated 
otherwise. Model parameters can be found in Table 1. 

(1) Resolution. In model sl6a, we test the convergence properties of our numerical 



^The ratio of vertical mass loss rate to the mass accretion rate is 



3.4. Dependence on Model Parameters 



Mr 




(16) 



where Mr is the mass accretion rate at the disk radius r, 6Y, is the change of disk surface density in 6t, and 
we have used a disk turbulent viscosity v — ac^/il. For r/H = 30, Mz/Mr ^ 1. 
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models by doubling the resolution of the fiducial run to 768 x 512 x 256. We run this model 
to ~ 75 orbits. Averaging over the last ~ 25 orbits in the satuated state, we found 
that (a) ~ 0.023, almost double of that in the fiducial run {{a) ~ 0.013); at the highest 
resolution explored in this work, saturation level continues to increase with resolution. The 
dynamic range in resolution explored in this work is modest (highest resolution in this work 
is ~ 20 — 40 zones per H) due to the computational demands of the large bojjfl- We have also 
monitor ed "quality f a ctor " Q, where Qi = AMRi,i/Axi = 27rt>a,i/(fiAxj) (see a discussion 



of Q in iNoble et al.l ( l2010l )). the zones per most unstable linear MRI wavelength in our 
calculations. For the bulk of the disk inside ±2H region, our highest resolution run gives a 
volume and time averaged Qy = 33.7 and = 6.8 in the saturated state, and so by this 
measure the toroidal field MRI is well-resolved while the vertical field MRI is marginally 
resolved. Of course, the evolution of the disk is not well described by linear theory in the 
fully turbulent state, so it is not clear whether Q is a good indicator of when MRI-driven 
turbulence is sufficiently resolved. 

It is worth mentioning the convergence properties of shearing box simulations done in 
smaller boxes: (a) the unstratified box with a net toroidal field, (b) the stratified box with 
periodic vertical boundary conditions and (c) the stratified box with outfiow boundaries. 
First, using asimilar algorithm in unstratified disk simulations with a mean toroidal field. 



Guan et al.l (l2009l ) reported that in the resolution range 32 — 256 /H saturation energy in- 

1 /3 

creases with resolution (oc Nx ). They also pointed out that convergence is expected at 
higher resolution when the energy containing eddies are resolved. For the stratified disks. 



(IShi. Krolik. &: Hirosdl2010l ) used L^. = 2H stratified shearing boxes simulations with ver- 
tical outfiow boundaries and they also found that (a) increases with resolution and an 
(a) ~ 0.035 in their highest resolution run at 32/ H. Recently, stratified disk simulations 
done in smaller boxes [Lx ~ H) and periodic boundary conditions with zero-net fiux have 
demonstrated convergence with (a) ~ .01 with a resolut ion ~ 32 — 128 /H using ATHENA 



code and periodic boundary conditions ( jPavis et al.ll2010l ). The sustained turbulence may 



be due to the presence of a mean toroidal field in the disk midplane. Notice that (a) in 
their work is normalized with the initial midplane pressure Pq, which is normally a factor 
of a few larger than domain-averaged (P) used here. Using the definition in this work, 
their(a) ~ 0.04. It is therefore possible that in stratified disk simulations, net-toroidal-field 
and zero-net-fiux models will have similar convergence properties. If this is the case, we then 
expect a convergence at 64 — 256 /H using our ZEUS-type cod^. 



■^For example, a (Lx, Ly, Lz) = (16, 20, 10)-ff run with resolution 768 x 512 x 256 (run sl6a) and tf ^ 150 
orbits required ^ 0.5 x 10^ cpu hours on abe cluster at NCSA. 

run of our fiducial run size with a resolution 6A/ H would require > 5 x fO^ cpu hours on NCSA's abe 
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(2) Lx- For we have carried out runs with = H,8H, and 32H, denoted by si, 
s8, and s32. Time averaging the last 50 orbits in each run, we found the saturation level in 
all these runs are close, with (a) ~ 0.0191 ± 0.00453 when = H^Ja) 0.0124 ± 0.00116 
when L.^ = 8i/, (a) ~ 0.0125 ± 0.000965 when = 16H, and (a) ~ 0.0269 ± 0.00211 when 
Lx = 32H, where the numbers after ± denotes standard deviation a. The dependence of (a) 
on the box size is not clear, however, it is difficult to measure (a) in the = H box because 
of the large fluctuations. In runs with > SH, the cr-to-mean ratio is aroud 0.07 — 0.09, 
while Lx = H gives a a-to-mean ratio ~ 0.25. 



Past stratified disk studies ( (jPavis et al.ll2010l : IShi. Krolik. &: Hirosell2010l )) have shown 
that there exist significant (order of unity) long term fiuctuations m. ^ H box. The 
evolution of magnetic energy density in the disk {EB,d) for = H and = IQH runs are 
shown in Figure [Hi The smaller fiuctuation in large models suggest that: (a) parts of the 
disk with horizontal separation > H are uncorrelated, and (b) the volume integration over 
large-enough domain will smooth out these local fiuctuations. 

What have we learned in these lar ge domain size rnodels with > lOiJ? Our = H 
run is similar to the toroidal model of iMiller fc Stond (120001 ). First, in large box runs, the 
plane-averaged vertical disk structures are similar to those in smaller box runs: we have 
observed a gas-pressure supported disk with a Gaussian density profile inside ~ 2H and an 
extended magnetic dominated corona outside ~ 2H. Second, the long-term average of disk 
turbulence saturation level is also very similar to the ~ H runs, albeit with much smaller 
temporal fiuctuations. Statistically, for saturation measurement purposes, a large domain 
run can be regarded as a sum of smaller H run, where the temporal and spatial fiuctuations 
are smoothed out by integrating over decorrelated disk regions. 

Our models also suggest that a magnetically dominated corona cannot be studied in 
an Lj, ~ if box (if it can be studied in a numerical MHD model at all). In large domain 
size runs with > lOH at \z\ > 2H we find features in the magnetic field correlation 
function on scales of ~ lOH, indicating the existence of meso-scale structure. Although 
the magnitude of the mass, angular momentum and energy transport in the corona is small 
compared to that in the central disk, the corona and the central disk are dynamically con- 
nected and large scale structure in magnetically dominated upper layers may still infiuence 
the spatial correlations/structures of the disk below (we will explore this issue in a forthcom- 
ing paper). Therefore, in accretion disk models where the spatial structure of the corona is 
importa nt, such as phenomenological models for accretion disk coronae (e.g., the statistical 
model of lUzdensky &: Goodmarul2008l ) and disk spectra calculations, the radial extent of the 



cluster. 
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corresponding numerical simulation may require a > IQH. 

(3) Lj,. We have investigated the effect of vertical boundaries by running a model (sl6c) 
with Lz = 12H. We find no qualitative difference between the = 12H and = lOH 
models: the saturation (a) ~ 0.0141, jE^/pocl ~ 0.0125, and {Eb,c)/pqcI ~ 0.00497, all 
within ~ 10% of the fiducial model. We also obtain a similar vertical disk structure when 
increases: inside 2.bH the disk has a well-fitted Gaussian [p\{z) and a flat [Eb\{z)] outside 
2.bH the corona extends to \a\ = ±6H in model sl6c. Both [p]{z) and [Eb]{z) also have an 
exponential profile, but the fitted coronal exponential scale height increases ~ 20% compared 
to the fiducial model. The coronal loop distribution functions are almost identical to those 
of the fiducial model, which is not suprising considering the steep decline of d^/dAzj^s.^ as 
^^max exceeds ~ H. As discussed before, some caution is needed in interpreting this vertical 
extension of corona structure with increasing domain size: our calculation is essentially a 
MHD calculation, whereas real disk coronae are probably force-free and also infiuenced by 
non-ideal plasma effects (e.g., reconnection) that are ill- modeled in our numerical scheme. 

(4) I3q. We have tested the effect of initial field strength on the saturation level. In most 
of our runs we start from a uniform toroidal field inside the disk with /3o = 25. We then carry 
out a comparison run sl6b with the same field geometry but weaker strength /3o = 100. We 
find that turbulence saturates at the similar level using weaker initial field strength, with 
Jo) ~ 0.0157, {Eb,a)/PocI ~ 0.0152, and {Eb,c) I PoC^ ~ 0.00647, therefore in stratified disks 
the saturation does not depend on the initial field strength. 

In comparison, in unstratified disk model s (Eb) is found to scale with the initial mean 
field strength {ByjQ (see a detail discussion in lGuan et al.l 120091 ). The important difference 



here lies in the stratification and the accompanying outfiow boundary conditions, which allow 
changes in mean toroidal field strength in the turbulent disk. The stratified disk model then 
allows the disk to adjust its net fiux and field strength in a self-consistent way. It is worth 
pointing out that the saturation level is a volume average over large scales where different 
parts of the disk decorrelate; on the scale where the turbulence is localized (< if), it is still 
possible that the local saturation (-EB,d)iocai (-^j/.d)iocar 

Does the saturation level depend on the instantaneous mean field strength in the strat- 
ified disk? The evolution of the mean azimuthal field (-Bj/,d) in the region \z\ < 2H in model 
s32 is plotted in Figure HI The mean field does not have a fixed value and changes signs over 
a time scale of ~ 10 orbits. Averaging the last 100 orbits, the mean magnitude of toroidal 



iGuan et al.l (|2009l ) found a linear relation between (By)^ and the saturation a oc (Eb) oc PqCsVa.j/Oi 
where VA,yo = Byo/li-Kpociy/^ is the mean azimuthal Alfv en speed. This result is also consistent with 
scalings obtained from earlier work fe.g.. iHawlev et all (jl995l) ). 
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field strength in the turbulent disk region is \ {By) \ ~ 0.012\/47rpQ^^Cs, therefore there is weak 
net toroidal field in the turbulent disk region to drive MRI. The figure also shows the evo- 
lution of the mean magnetic energy density {Eb,a} and mean xy stress {Wxy,d) in the same 
run. Both {Eb^a) and (Wxy,<i) have a fixed overall saturation level with a superimposed small 
oscillation with a period half of that of {By). Again, this is dramatically different from the 
unstratified disk, where the saturation level is proportional to the mean field strength. The 
overall saturation level in a stratified disk is not determined by the instantaneous mean field 
strength, nor by the initial field strength. 

On the other hand, as shown in Figure [HI the oscillations of {By) and {Eb^a) are closely 
correlated and oscillation period for {By) is twice that of {Es^d)- This suggests that {B'^) is 
correlated with (-By), even though {By) <^ (5^)^^^. Therefore, the saturation level may be 
determined both by the MRI induce d by the mean to roidal field in the turbulent disk region 



and magnetic buoyancy effects (e.g.. IVishniad (120091 )). 



4. Butterfly Diagram: A Mean Dynamo In the Disk? 

One interesting feature appearing in all our models is an oscillation of the mean magnetic 
energy on a timescale of a few orbits. As an example, we plot the "butterfly" diagram for 
model s32 in Figure [T0| which illustrates the evolution of [Eb\{z). This bears a superficial 
resemblance to the famous butterfly diagram observed in solar activity cycles. 

We use Fourier analysis to determine the period of butterfly diagram. Using data 
j dyEB{\z\ = 2.5H] x,t), taken from two layers with \z\ ~ 2.5H and have been averaged in 
y direction to improve statistics, we perform a two dimensional FFT (in x and t) on the 
data set. The normalized temporal power spectral density (PSD) for [Eb^\z\=2.5h] in model 
s32 are shown in Figure [TTJ Here we have plotted a cut through = plane in the 2D 
kx — f PSD map. We have also checked that the different sides of the disk have very similar 
PSD and we have plotted the sum of contribution from both layers. The arrow in the figure 
marks the peak frequency in the PSD. This frequency, / ~ 0.03fi corresponds to the period 
of the butterfly diagram for [Eb]- 

The PS D has P ^ f^, with k ~ —2.3. Interestingly, results from recent global GRMHD 



simulations (iNoble fc Krolikll2009l ) suggest that the slope for the coronal luminosity temporal 
power spectrum is ~ —2, almost independent of model parameters and very close to what 
has been observed at high frequency in black hole accretion disk systems. The power- 
law index for the temporal power spectrum from local and global simulations are therefore 
remarkably close, considering we are only calculating the temporal spectrum for coronal 
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magnetic energy density. 

The period for [Eb] is Pieb] ~ 5 orbits. Besides [Eb], one could also plot butterfly 
diagrams for [Wxy] and [By]. The period for [W^y] is the same as that of [Eb], ~ 5 orbits, 
while the period for [By] is twice that of [Eb], P[By] ~ 10 orbits, because of the reversal of 
mean flelds (see Figure [9]). 

For [Eb], we flnd P{Eb) ~ 5 orbits in all our models. This quasi periodicity has appeared 
in all the stratifled shearing box simulations that we are aware of, even in those with p e riodic 



vertical boundary conditions (Stone et al., 2009). Interestingly, [Reynolds fc FabianI (|2008[ ) 
has also obtained similar butterfly diagrams at certain radii (e.g., r = Sr^ and r = lOr^, 
where = GMj(? is the gravitational radii) in their global pseudo-Newtonian thin disk 
simulations. It would be of interest in the future to test (a) whether the butterfly diagram 
is simply a local feature at a certain location on the disk (as in shearing box simulations) or 
this quasi-periodicity can be coherent and sustained over a large radial range, and (b) what 
model parameter(s) the period depends on. 

The butterfly diagram together with the reversal of the mean fields (for both the dom- 
inant toroidal field and a weak r adial field) in the disk may be modeled by a mean field 
dynamo of a type (e.g., iMoffatt (|l978|)@). In the rest of this section, we will present a toy 



model to give a qualitative description of these oscillations. 

Let us first consider two important dynamical processes in a stratified disk: (1) the MRI- 
driven turbulence, which draws free energy of rotation and operates on the orbital timescale 
~ fi^^; (2) magnetic buoyancy, which operates on the local Alfven timescale ta ~ 
In our simulations we found in the disk region \z\< 2H the magnetic energy density is almost 
a constant with height, with [Eb] > lO^^poC^, which gives an average magnetic buoyancy 
timescale ta < a few Vt^^ inside the disk. The period of the butterfiy diagram is much 
longer than these two timescales. Therefore these two processes alone can not describe the 
dynamics represented in butterfiy diagrams. 

The a type mean dynamo equations for the disk mean fields (-B^.) and {By) can be 
derived from averaging the induction equation, dt{B) = (V x (i? x B)), here ( ) denotes 
ensemble averages. Assuming the turbulent EMF e is related to the mean field with a 
dynamo parameter (e) = {6v x 6B) = ai{Bi), one simple form of dynamo equations in 



^In this work we use a to denote dynamo model type. It should not be confused with the accretion disk 
turbulence level parameter a 
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a stratified thin Keplerian disk is (cf., Eqn (5-6) in IVishniacfc Brandenburgl ( 119971 )) 



^tW = -^^(5x)-5.(K)W) + 9,(ai(5,)) , (17) 

and 

dtB, = -d,i{v^){B^)) - d^ia^iBy)) , (18) 

where is a characteristic vertical velocity induced by magnetic buoyancy. In Eqn f[T7|) the 
first term is the shear term, second term denotes buoyancy due to the mean field, and the 
last term is the mean field dynamo term. Only dz terms are retained because the disk is 
thin. For simplicity we have also dropped the diffusion terms. We then take dz ~ 1/{2H) 
and fb ~ \va\ = \B\/^/4^Fpo ~ iByl/y/Anpo, where \va\ is the mean Alfven speed. Eqn lfTTl) 
and Eqn fll8p then become 

^ = --nB,^-^-^By + ^B^, (19) 
dt 2 2H ^ 2H ' ^ ' 

and 

dB^ \va I „ tt2 „ /X 

— ^ = -^B^ -B^ . 20 

dt 2H 2H y ^ ' 

For clarity we have dropped ( ) in the above equations. Notice that Eqn (JT9l) and Eqn (1201) 

have no spatial dependence. Taken together, they are coupled ODEs and can be solved 

numerically given initial conditions for By and B^. 

In Figure [I2] we plot one solution for this toy model. This solution is obtained by 
integrating the above equations from an initially pure toroidal field with /3o ~ 22 and by 
choosing ai = 6:2 = —0.010. The period for By in this particular toy model is ~ 10 orbits. 
The magnitude of a controls the oscillation frequency: in general, larger |a| leads to smaller 
period, although the scaling is not linear. Initial conditions have little effect on the evolution 
in our toy model. In conclusion, the butterfly diagram and the mean field reversal observed 
in these simulations may imply a mean field dynamo at work in stratified disks. 

Does it make sense to identify these oscillations with observed QPOs? In a Keplerian 
disk the orbital frequency at r is /orb = l/(27r)(G'M)^/^r~^/^. The QPO frequency is /qpo = 
1/5/orb ~ 20x (r/10M)"^/^(M/10MQ)"^Hz. Our disk model represents a geometrically thin, 
optically thick disk. This is most easily understood as corresponding to the high soft state in 
black hole X-ray binaries, which is dominated by a thermal component. For a 10 black 
hole a 5Hz QPO (e.g., XTEJ1550-564) corresponds to r^^ ~ 25M, which is far from innermost 



''By definition, ai — — (^^'v^^^^ SvjBy) ^ _ _ {Sv;,sb^^ Sv^sb^) ^ principle, ai does not 
necessarily equal a2 due to anisotropy. 
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region of a thin disk where most of the thermal X-ray emissions presumably originates. This 
oscillation frequency may be sensitive to the disk vertical structure (e.g. if the disk is not 
isothermal), and therefore may exhibit a much more complex behavior in real disks, in 
which the vertical structure is closely coupled to vertical energy transport. On the other 
hand, observations indicates that QPOs are absent or very weak in the thermal state, but may 
appear in the very high state when a sizable thermal disk c omponent is present, although t he 
QPOs are more associated with Comptonizing electrons (( iRemillard fc McClintockl 120061 )): 
it is difficult to associate the butterfly oscillations with observed QPO phenomena. 



5. Summary and Discussion 

We have carried out stratified shearing box simulations with domain size = H to 
Lx = 32H to study properties of isothermal accretion disks on a scale larger than the disk 
scale height H. Our numerical models have vertical extent > 5H above and below the disk 
midplane with outfiow boundary conditions. All models start from a net mean toroidal field 
in the central disk region and the mean fields are allowed to change in the evolution. 

We find the disk has an oscillating mean toroidal field and (a) ~ 0.012 — 0.025 in the 
parameter range we explored. We have not found a clear dependence of (a) on in our 
models, although the temporal variances in volume averaged quantities decreases with L^. 
The highest resolution used here is modest (20 — 40 zones per H), and we have observed (a) 
increases with resolution. Recently, Stone el al. report a converged (a) ~ 0.04 in L^. ~ if 
high resolution stratified disk simulations with zero-net-fiux and periodic vertical boundary 
conditions (so that the volume-averaged field cannot change during the evolution). The 
sustained turbulence may be due to the presence of a mean toroidal field in the region close 
to the disk midplane, lending plausibility to the idea that the saturation mechanism of MRI 
in stratified disks near the midplane is similar to that in unstratified disks with a net toroidal 
field. 

In the saturated state the disk vertical structure consists of (a) a turbulent disk at 
\z\ < 2H and (b) a magnetically dominated upper region at \z\ > 2H, confirming earlier 
small {Lx ~ H) box results. 

At \z\ < 2H, the disk is mainly supported by gas pressure, and a Gaussian density 
profile is observed. The plane averaged magnetic energy density [Eb] (z) and Maxwell stress 
[Mj.y](z) are nearly uniform with vertical height z in this region, where the disk is marginally 
stable to the Parker instability. At \z\ > 2H, exponential dependences on z are observed for 
both [p] and [Eb]. Fitting formulae for [p]{z) and [Eb]{z) are given in Eqn and Eqn (IHl) 
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respectively. 

Using a two-point correlation function analysis, we found that close to the midplane, the 
disk is dominated by small scale (< H) turbulence, very similar to what we have observed 
in unstratified disk models. In the corona, magnetic fields are correlated on scales of ~ lOH, 



implying the existence of meso-scale structures. Recently iJohansen et al.l (120091 ) have also 
observed large scale pressure and zonal flow structures in their large shearing box simulations. 
We will give a detailed report of meso-scale structure in isothermal disks in a forthcoming 
paper. 

We have adopted a statistical approach to study the geometry of coronal magnetic fields. 
Only ^ 4% of coronal field lines are open. For closed field lines, we calculated the magnetic 
loop distribution function for the loop foot separation Ax in the x — y plane, loop maximum 
height A2;jnax5 and loop orientation angle 6'foot- The loops are dominantly toroidal due to 
the differential shear. The loop foot distribution between H — 20H is a power law with 
an index /c ~ — 1.25. In the phenomenological model of UG, this corresponds to the limit 
where reconnection is slow compared to the shear. These comparisons are limited because 
our models are working in an ideal MHD regime and reconnection is purely numerical. 

In our models both vertical energy and momentum flux are negligible in the steady 
state. The mass loss rate from the disk surface is small and decreases with increasing L^. 
The surface effects are therefore minimal and indicate a lack of disk winds in our stratified 
disk models. The weak winds are consistent with the constraint that we have a zero-net 
vertical magnetic flux in these models. A Blandford-Payn e type wind requires the existence 



of a vertical net field (e.g., see ISuzuki fc Inutsukal (120091 )). although we note that in their 



models the most unstable wavelength for the extremely weak field are probably not resolved.) 
Initial investigations show that even a weak (/3o ~ 1600) net z field will induce very violent 
accretion in stratified shearing box models: at certain region of the disk accretion will run 
away, eventually causi ng the disk break into rings. Similar phenomena were reported in net 



vertical field models of iMiller &: Stond ( 120001 ). 



We have confirmed the "butterfly" diagram seen in earlier stratified disk models of size 
~ H. The butterfly diagrams persist even in our largest runs with = 32H. We 
also report the reversal of the mean fields (for both the dominant toroidal field and a weak 
radial field) in the disk on a timescale twice that of [Eb]- The periods for the butterfly 
diagram are close in all our models, P ~ 5 orbits for [Eb] and P ~ 10 orbits for [By]. 
The mean field reversal and butterfly diagram may indicate the existence of a mean field 
dynamo in stratified disks, perhaps controlled by the MRI and magnetic buoyancy. We 
have presented a toy model for an a type mean field dynamo in stratified disks and found an 
ttimp ~ 0.01 will produce the reported period. Further exploration of parameter dependences 
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would be useful for analytical modeling. In the future it would also be interesting to test 
whether the butterfly oscillations persist when averaging over a large range of radii in global 
disk simulations. The butterfly diagram may be associated with low frequency QPOs and 
therefore a good observational diagnostic for accretion flows. On the other hand, we also 
report a power-law index k ~ —2.3 in the temporal power spectrum for coronal magnetic 
energy fluctuations, consistent with results from recent GRMHD black hole accretion disk 
simulations. 

Our stratified disk models are primarily limited by the assumption that the disk is 
isothermal. Effects of thermodynamics and radiation therefore are neglected in this work. 
Our models are also limited by finite resolution, box size, evolution time, and the absense 
of explicit dissipation. Additional insights may also provided by the future explorations on 
magnetic field strength and geometry in disks. 
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Table 1. Model Parameters 



Model 


Size 


Resolution 




(a) 


{Eb,\z\<2h}/POCs 


{EB,\z\>2H}/p0c'i 


stdl6 


(16,20, m)H 


384 X 256 X 128 


25 


0.0125 


0.0121 


0.00427 


sl6b 


(16,20, 10)H 


384 X 256 X 128 


100 


0.0157 


0.0152 


0.00647 


sl6c 


(16,20,12)// 


384 X 256 X 160 


25 


0.0141 


0.0125 


0.00497 


si 


(1,20, 10)H 


48 X 256 X 128 


25 


0.0191 


0.0171 


0.00933 


s8 


(8,20,10)// 


192 X 256 X 128 


25 


0.0124 


0.0115 


0.00665 


s32 


(32,20, 10)// 


768 X 256 X 128 


25 


0.0269 


0.0270 


0.0106 


sl6a 


(16,20, 10)// 


768 X 512 X 256 


25 


0.0230 


0.0181 


0.0101 
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Fig. 1. — Snapshots of density and magnetic energy density in the fiducial model, taken at 
t = 100 orbits. Left: density p; right: magnetic energy density Eb] top: image at ?/ = 
plane; middle: image at x = plane; bottom: image at 2; = plane. 
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Fig. 2. — Vertical profiles of several x — y plane averaged quantities in the fiducial model. 
Upper left: density; Upper right: magnetic energy density; Lower left: plasma (3; Lower 
right: xy component of Maxwell stress M^y = —BxBy/An (solid lines) and Reynolds stress 
Rxy = pVxSvy (dotted lines). All quantities in solid and dotted lines are averaged from the 
last 50 orbits. To illustrate the time average effect, we also plot [Eb]{z) at t = 900^2^^ 
(dashed lines) in the upper right panel. The slight asymmetry of [Eb]{z) and [M2.,y](z) in 
the 1^1 < 2i? region is probably due to our choice of orbits interval for time average. 
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Fig. 3. — Evolution of (a) (solid lines), (Eb^) / Pocj. at \z\ < 2H (dotted lines), and 
{Eb,c)/ Pocl at \z\ > 2H (dashed lines) in the fiducial model. Saturation (a) ~ 0.013 when 
averaged over the last 50 orbits. 
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Fig. 4. — Contour plots of 2D two-point correlation function ^b{z) for SB. Plotted are 
^B{z)/{4:7rpocl) in the (Ax, Ay) plane at three different vertical heights in the fiducial model. 
Left: mid plane; Middle: z = 2.5H; Right: z = 4.5H. The contours run linearly from 
—0.058 to 0.229 for 20 levels; sohd lines: > 0; dash lines: < 0; the heavy line is the 
contour. 
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Fig. 5. — Magnetic field lines originating from the plane z = 2.5H in the fiducial model 
at t = 600Q~^. The lines are evenly sampled spatially from the x — y plane. Both the line 
width and the color (see the online version for a color version of this plot) denote the flux 
carried by each line at the footpoint, normalized by the total flux from the z = 2.5H plane. 
Majority of the field lines return to 2; = 2.5H, forming closed loops. 
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Fig. 6. — Loop angle distribution functions in the fiducial model: 6 denotes the angle 
between the foot separation vector and the y axis. Loops are stretched in the azimuthal 
direction due to the shear. 
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Fig. 7. — Loop-distribution functions in the fiducial model. Left panel: foot point separation 
distribution. Left curve is for Ax and Right curve is for Ay. The distribution for AL = 
(Ax^ + Ay^)^/^ almost overlaps with the Ay curve. The heavy line indicates a A; = —5/4 
slope. Right panel: loop height distribution. 
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Fig. 8. — Evolution of the disk magnetic energy density {EB,d} in the model stdl6 (solid 
lines) and si (dotted lines). Plotted are the mean magnetic energy density in the region 
1^1 < 2//. 
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Fig. 9. — Evolution of {By) (solid lines), {Eb) (dotted lines) and (yV^y) (dashed lines) at 
\z\ < 2H in model s32. The oscillation period for (By) is twice that of {Eb) and {Wxy). 
These temporal oscillations may be caused by a mean field dynamo. 
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Fig. 10. — Butterfly diagram for [Eb] in tlie model s32. 
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Fig. 11. — Normalized temporal power spectral density for [Eb] in the model s32. The data 
are taken from the layers with \z\ ~ 2.5H. We also draw the best-fit k — —2.3 slope for the 
temporal PSD. The arrow marks the peak frequency in the power spectrum. This frequency, 
/ ~ 0.03O, corresponds to the period of the butterfly diagram for [Eb]- 



-36- 




Fig. 12. — Evolution of By (solid lines) and (dotted lines) in our mean field dynamo toy 
model, ai — a2 — —0.01 in the plotted model. 



